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I.  INTRODUCTION 


This  report  describes  a  series  of  numerical  experiments  based  on 

a  computer  program  developed  at  the  General  Electric  Company  by  Scala 
1* 

and  Gordon  .  This  program  which  is  based  on  the  full  time -dependent 
Navier-Stokes  equations  is  used  to  calculate  the  flow  field  around  a 
circular  cylinder  accelerated  to  a  supersonic  velocity.  The  validity 
of  our  numerical  results  is  determined  by  comparison  with  analytical 
solutions  and  experimental  wake  measurements. 

Accurate  numerical  solutions  of  the  full  Navier-Stokes  equations 
are  desirable  for  many  problems  in  weapons  technology.  For  problems 
like  the  internal  flow  in  fluidic  devices  or  the  wake  flow  behind  a 
ballistic  missile  the  simplifying  assumptions  of  inviscid  flow  or 
high  Reynolds  number  'boundary -layer  flow  are  not  applicable.  To 
describe  the  flow  field  in  such  problems  it  is  possible  to  resort  to 
numerical  solution  of  the  full  Navier-Stokes  equations. 

Several  numerical  solutions  of  the  Navier-Stokes  equations  iiave 

been  published  for  compressible  viscous  flows,  such  as  those  of 

2  3  4  5 

Crocco  ,  Thommen  ,  Scala  and  Gordon  ,  Allen  and  Cheng  ,  and  Trulio 

and  Walitt^.  The  flow  detail  exhibited  in  these  solutions  is  quite 

remarkable,  but  in  certain  regions  of  the  flow  field  the  numerical 

solutions  may  not  be  very  accurate.  At  the  present  time  we  should 

accept  the  results  of  any  such  numerical  computation  with  great 

caution  since  the  mathematical  theory  on  which  the  numerical  solution 

is  based  is  inadequate.  There  are  no  rigorous  proofs  of  stability  or 

convergence  of  the  numerical  methods  employed  for  these  solutions,  nor 

are  there  any  proofs  of  uniqueness.  It  is  necessary  to  assess  the 

validity  of  any  such  numerical  computation  by  comparing  it  against 

known  analytical  solutions  and  experimental  measurements . 

Results  are  presented  for  Mach  two  flow  past  circular  cylinders 
at  Reynolds  numbers  1J.0,  46.8,  157*2,  and  70A.6.  The  numerical 
solutions  are  compared  with  analytical  solutions  for  the  density  ratio 


*  M.e  tinted  on  page  ^ 
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across  the  how  shock,  shock  structure,  and  shock  standoff  distance. 

The  validity  of  the  cylinder  wake -region  calculations  is  determined  by 

8 

comparison  with  experimental  hot-wire  wake  measurements  . 


II.  FORMULATION  OF  THE  PROBLEM 

The  governing  equations  on  which  our  numerical  solutions  are  based 
are  formulated  within  the  framework  of  classical  continuum  mechanics. 
The  Navier-Stokes  equations  are  formed  from  Cauchy's  First  Law  of 
Continuum  Mechanics 


Dt 


3t 


=  g 


dx , 


(1) 


when  the  stress  tensor  t.  is  expressed  as  a  linear  function  of  the 

i-0 


deformation  gradient.  We  thus  assume 

du  du. 

'  -  <  -  * +  x  sr-> 


du . 


dXj 


-)  • 


(2) 


The  basic  equations  have  been  written  here  in  Cartesian  tensor 
notation  for  simplicity.  The  viscosity  coefficients  appearing  in 
equation  (2)  are  further  assumed  to  be  related  by  Stokes  postulate 

3X  +  2p  =  0  .  (5) 

Equation  (3)  can  be  shown  to  be  valid  for  monotomic  ideal  gases. 

The  Navier-Stokes  equations  give  us  three  equations  to  relate  the  six 
field  functions.  We  also  use  the  continuity  ecuation 

au. 


52.  +  p  =  0  , 

Dt  H  dx.  * 

l 

the  energy  equation 

UP  ,  „  Sui 

pcv  Dt  +  P  dxT 
1 


(4) 


d_ 

dx 


/,  dT  \  .  /Su1n2 

^  dx.  dx, 

i  x  l 


du .  du.  du 

,  i  /  i  ,  _ m 

+  ^  S3T  '  5x~  +  5x7 

m  m  i 


00 


)  , 


1U 


and  the  equation  of  state  for  a  perfect  gas 

p  =  pRT  .  (6) 

In  this  study  the  time -dependence  in  equations  (l),  (U),  and 

(5)  is  retained  and  we  seek  steady-state  solutions  of  the  governing 
equations  by  time -wise  integration  from  a  prescribed  initial  state. 

The  steady  state  is  thus  approached  asymptotically  with  time. 

The  partial  derivatives  in  the  governing  equations  (l)  thru 

(6)  are  replaced  by  differences  in  both  space  and  time.  This 
approximation  reduces  the  system  of  nonlinear  partial  differential 
equations  to  a  system  of  nonlinear  algebraic  equations  which  can  be 
solved  by  numerical  iteration.  The  set  of  difference  equations  used 
in  this  procedure  was  developed  by  Scala  and  Gordon"'’. 

The  complete  system  of  equations  is  formulated  in  general  two- 
dimensional  orthogonal  coordinates.  This  system  of  equations  together 
with  their  difference  approximations  are  listed  by  Scala  and  Gordon1. 
In  the  present  study  we  consider  time -dependent  flow  past  a  cylinder 
using  cylindrical  coordinates  (R,9)  with  0  measured  from  the  forward 
stagnation  point. 


III.  DIFFERENCE  EQUATIONS 

The  complete  set  of  difference  equations  used  in  v-he  computer 
program  to  approximate  the  governing  equations  are  discussed  in 
detail  by  Scala  and  Gordon1.  The  procedure  followed  is  to  split  the 
full  system  of  governing  equations  into  "hyperbolic"  subsystems  and 
parabolic  terms  which  are  differenced  separately  and  then  combined 
into  a  final  set  of  difference  equations. 

The  second-order  viscous  terms ,  the  parabolic  terms,  are- 

9 

differenced  using  an  alternating  explicit-implicit  scheme  .  At 
every  time  step  the  calculation  at  each  interior  point  is  alternated 
between  an  implicit  formula  and  an  explicit  formula.  According  to 
this  procedure  implicit  equations  are  used  for  the  calculations  at 
every  alternate  interior  mesh  point  (with  respect  to  both  R  and  0)  at 
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time  t  =  (n)  At.  At  the  remaining  interior  mesh  points  explicit 
equations  are  used  for  time  t  =  (n)  At.  A  time  t  =  (n+l)  At  the 
explicit  points  of  the  previous  time  step  become  implicit  and  vice 
versa.  This  procedure  then  repeats  for  subsequent  time  steps. 

The  "hyperbolic"  subsystems  are  made  up  of  both  the  inviscid 
terms  and  the  lower  order  viscous  terms.  These  "hyperbolic"  terms  are 
differenced  by  transforming  them  to  diagonal  form  and  then  differencing 
them  according  to  the  sign  of  the  characteristics.  The  calculations 
for  these  terms  are  then  alternated  between  explicit  and  implicit 
equations  in  the  same  manner  a?  described  for  the  second-order  terms. 

The  set  of  difference  equations  used  in  these  calculations  and 
the  manner  in  which  they  are  numerically  integrated  introduced  several 
uncertainties  into  the  accuracy  of  the  numerical  solutions.  In  the 
present  difference  formulation  the  first-order  nonlinear  terms 
involving  spatial  derivatives  of  h(t)  k(T)  in  the  momentum 
equations  are  treated  like  hyperbolic  terms.  These  terms  do  not 
necessarily  remain  hyperbolic  throughout  the  flow  field.  The  amount 
of  error  inbroduced  into  the  steady-state  solutions  by  this  approximation 
cannot  be  evaluated. 

The  second  uncertainty  arises  from  the  implicit  evaluation  of  the 
second-order  cross-product  derivatives  (d^/dRde  ,  etc . ) .  These 
derivatives  are  approximated  at  alternate  time  steps  by  assuming  that 
implicit  quantities  are  known  from  the  previous  time  step.  This 
simplification  eliminates  the  iterations  needed  if  the  true  implicit 
equations  were  used  at  es.ch  alternate  time  step.  This  procedure 
introduces  some  uncertainty  into  the  time -dependence  exhibited  by  the 
numerical  solution.  It  is  not  known  what  effect,  if  any,  this  procedure 
has  on  the  final  steady-state  numerical  solution  which  is  obtained 
asymptotically  with  time. 
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IV.  INITIAL  CONDITIONS  AND  BOUNDARY  CONDITIONS 


The  treatment  of  the  time -dependent  Navier -Stokes  equations 
requires  the  specification  of  both  initial  conditions  (t=0)  throughout 
the  flow  field  and  boundary  conditions  (t  >  0)  on  all  boundaries. 

Figure  1  gives  a  schematic  representation  of  the  flow  field  surrounding 
a  cylinder  traveling  at  supersonic  speed.  The  flow  field  is  character¬ 
ized  by  a  region  far  upstream  of  the  body  in  which  viscous  effects  are 
negligible.  There  are  two  shock  regions  in  the  flow  field;  the  bow 
shock  and  the  wake  shock  regions .  Near  the  surface  of  the  cylind  »r 
viscous  effects  become  predominate  and  this  region  is  characterized  by 
a  viscous  layer  which  merges  into  a  wake  region  downstream.  In  the 
wake  region  neither  inviscid  nor  boundary-layer  considerations  are 
valid . 

For  times  s  0  in  our  numerical  computations  we  assumed  a  uniform 
flow  having  freestream  conditions  throughout  the  flow  field.  For  times 
t  >  0  the  velocity  components  on  the  surface  of  the  cylinder  were 
brought  linearly  to  zero  and  held  at  this  value.  The  following 
velocity  boundary  conditions  were  specified  on  the  cylinder's  surface 
u(0i)  =  Bua(ei)  ,  (7) 

v(ei)  =  evjep  ; 

(fcl  -  t)/t1  for  0  S  t  i  t.  , 

1  (8) 

0  for  t  >  t  • 

For  times  t  i  we  thus  specify  the  no-slip  condition  on  the  surface 
of  the  cylinder. 

Physically,  the  process  of  initialization  described  by  (7)  and 
(8)  treats  two  different  types  of  initial  conditions.  If  t^  is  less 
than  the  time  step  At  used  in  the  numerical  integration,  then  we  model 
a  cylinder  impulsively  accelerated  from  rest.  If  t  >  At  we  treat 
flow  past  a  porous  cylinder  with  the  amount  of  mass  flux  through  the 
surface  decreasing  to  zero  between  times  t=0  and  t=t^. 


where 
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•rounding  Circular  Cylinder 


In  the  latter  ease,  for  times  t  <  t^,  mass  exits  the  downstream 
surface  of  the  cylinder  into  the  flow  field.  We  specify  the  density 
of  this  mass  to  be 

P=|p(  3+  (3)  (9) 

1 

for  times  0  s  t  <  t^  .  The  procedure  followed  by  Scala  and  Gordon 
for  specifying  density  boundary  conditions  is  that  whenever  the  flow 
enters  the  computation  region  through  a  boundary  the  density  must  be 
specified  a3  a  boundary  condition,  otherwise  it  is  calculated  at  that 
boundary.  The  results  of  our  calculations  show  that  this  boundary 
condition  may  be  incorrect  but  it  is  not  known  whether  the  problem 
is  properly  posed  as  stated. 

The  boundary  condition  for  the  temperature  on  the  cylinder' s 
surface  is  specified  by  the  assumption  of  an  adiabatic  wall 


It  is  hoped  that  the  time -dependent  history  of  the  solution 
resulting  from  these  initialization  procedures  will  have  a  physical 
counterpart.  It  is  assumed,  however,  that  the  resulting  steady-state 
solution  will  be  unique,  independent  of  the  Initialization  process 
prescribed.  This  assumption,  unfortunately,  cannot  be  verified 
analytically. 

In  our  computations  we  choose  the  steady  state  from  observations 
of  the  time -dependent  flow  in  the  bow  shock  and  wake  regions.  At 
time  zero  the  bow  shock  forms  adjacent  to  the  cylinder  and  then 
propagates  upstream,  gradually  reaching  a  standoff  distance  which 
does  not  change  appreciably  with  time.  Tlie  wake  flow  then  establishes 
itself  after  the  bow  shock  has  reached  its  steady-state  position. 

After  the  bow  shock  is  stationary  we  monitor  the  wake  flow  until  this 
flow  is  invariant,  or  almost  so,  after  several  hundred  time  steps. 

The  flow  field  at  this  stage  is  then  said  to  represent  the  steady-state 
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nolutlon. 

The  specification  of  boundary  conditions  far  from  the  cylinder 
poses  a  difficult  numerical  problem.  We  have  used  two  different 
techniques  in  this  study  in  an  attempt  to  determine  the  influence 
of  these  boundary  conditions  on  the  accuracy  of  the  numerical 
solution.  The  first  method  follows  that  used  by  Scala  and  Gordon"*-. 

We  specify  "artificial"  boundary  conditions  along  an  outer  boundary 
located  at  a  finite  distance  from  the  cylinder's  surface.  The  second 
technique  has  been  used  by  Sills  .  In  this  method  a  coordinate 
transformation  is  used  to  map  the  interval  Co,®]  of  R  onto  a  finite 
interval  [o  l].  The  proper  boundary  conditions  at  infinity  can  then 
be  specified  at  the  boundary  of  the  transformed  region.  Table  AI 
shows  which  method  was  used  for  each  calculation. 

In  the  former  case,  with  the  outer  boundary  located  at  a  finite 
distance  from  the  cylinder,  we  prescribe  the  outer  boundary  to  be 

=  R2  [l  +  1  (m-l)  xJ  -  5  (m-l)  x*4'  ]  ;  (ll) 


where 


x 


0 

rr 


(12) 


The  value  m  determines  the  degree  of  difference  mesh  stretching 
downstream  of  the  cylinder.  Our  calculations  were  carried  out  for 
m  =  2.5 .  With  this  value  of  m  the  outer  boundary  downstream  of  the 
cylinder  (6  =  l80°)  is  two  and  one -half  times  farther  from  the 
cylinder  than  at  0  =  0° .  Figure  2  shows  a  schematic  of  the  region 
of  computation  with  m  =  2.5*  Transformation  (ll)  facilitates  the 
use  of  a  fine  difference  mesh  upstream  in  the  bow  shock  region  and  a 
coarse  mesh  downstream  in  the  wake  region. 

For  those  calculations  in  which  the  mapping  was  employed,  the 
radial  coordinate  was  transformed  by 


5 


1  -  e  '  C'(R-Ro) 


(15) 


The  arbitrary  constant  a  is  a  stretching  factor  which  determines  the 
location  of  the  image  points.  This  transformation  permits  the 
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specification  of  boundary  conditions  at  infinity  while  allowing  the 
interior  grid  points  to  be  concentrated  near  the  cylinder.  It  should 
be  noted  that  the  use  of  this  transformation  results  in  a  solution 
which  is  inherently  asymptotic  a.t  the  outer  boundary.  As  long  as  the 
derivative  d/3§  remains  bounded  as  §  1  in  the  computational  plane, 

then  d/SR  -*  0  as  R  -»  ®.  This  inherent  asymptotic  behavior  may 
unnecessarily  result  in  an  overspecification  of  boundary  conditions 
in  the  numerical  computation. 

In  our  numerical  solution  we  have  made  a  symmetry  assumption  about 
the  flow  field  in  order  to  reduce  the  computer  storage  requirements  and 
computation  time.  We  assume  a  plane  of  reflective  symmetry  bisecting 
the  cylinder's  cross  section  and  numerically  treat  flow  in  the  upper  half 
plane  0  <  0  ^  tt  .  .  We  thus  specify  symmetry  boundary  conditions  along 
boundaries  at  0  -  0  and  0  *  l80° .  The  symmetry  conditions  prescribed  on 
these  boundaries  are 


||  =  0,  ||  =  0  and  V  =  0  .  (14) 

No  boundary  condition  was  prescribed  for  the  density  on  these  bound¬ 
aries  for  the  calculations  with  a  finite  outer  boundary  (ll).  This 
follows  the  procedure  used  by  Scala  and  Gordon^-.  For  the  calculations 
using  the  mapping  (13)  we  additionally  specified 


(15) 


along  the  boundaries  at  0  =  0°  and  6  -  180° . 

In  all  calculations  we  prescribed  freestream  velocity  components 
and  freestream  temperature  along  the  outer  boundary,  whether  using 
(11)  or  (13).  Free stream  values  were  specified  for  density  at  the 
outer  boundary  whenever  the  flow  at  a  boundary  point  was  directed  into 
the  computation  region  when  using  (ll).  For  the  calculation  using 
the  mapping  (ljj)  freestream  values  for  density  were  always  specified 


at  infinity. 
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Figure  3  shows  the  domain  of  computation  and  boundary  conditions 
for  two  of  the  cases  calculated  using  the  finite  outer  boundary  (ll). 

The  boundary  conditions  specified  on  each  boundary  are  shown  in  thin 
f igurc .  The  two  oute-*-  -„ndaries  shown  in  this  figure  correspond  to 
their  locations  for  t  .lculations  performed  for  Re^  =  46.8. 

V.  DESCRIPTION  OF  COMPUTATIONS 

Numerical  computations  were  performed  for  Mach  two  flow  past 

circular  cylinders  at  freestream  Reynolds  numbers  of  15*0,  1*6.8, 

157.2  and  704.6,  based  on  cylinder  diameter.  Appendix  A  contains  the 

dimensional  values  assigned  to  the  flow  field  variables  for  each  case. 

Each  case  will  be  referenced  by  its  approximate  freestream  Reynolds 

number.  The  four  cases  for  a  Reynolds  number  of  46.8  are  differentiated 

by  a  suffix  (i.e.  47-1,  47-2,  47-3  and  47-4). 

The  largest  cylinder  considered  had  a  diameter  of  approximately 

one-half  inch  (Re^  =  704.6).  The  calculations  for  this  case  represent 
d 

an  upper  limit  of  the  type  of  calculation  that  can  be  performed  in  a 
reasonable  amount  of  computer  time  using  this  computer  program.  All 
calculations  were  performed  on  the  BRIESC  II  computer  at  the  Ballistic 
Research  Laboratories. 

The  difference  mesh  was  constructed  so  that  the  angular  spacing 
remained  constant  at 

&Q  =  (  rr/30)  =  (l6) 

for  all  computations.  The  length  of  a  mesh  element  in  the  radial 

direction  varied  with  the  different  computations.  The  smallest  mesh 

-3 

was  employed  in  case  47-1  for  flow  past  a  3.056x10  ft  diameter 
cylinder  with  £R  =  1.050x10  ft  along  the  forward  stagnation  stream¬ 
line.  The  largest  mesh  occurred  in  calculations  such  as  case  47-3 
which  used  transformation  (13).  In  such  cases  the  mesh  elements 
farthest  from  the  cylinder  stretched  from  some  finite  radial  coordinate 
to  infinity.  The  radial  mesh  variation  was  chosen  by  specifying  a  in 
transformation  (13)  so  that  there  would  be  several  mesh  points  per 
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freestream  mean  free  path  across  the  bow  shock  at  all  times. . 

A.  Re^  =  704.6 

The  calculations  for  Rej  =  704.6  were  initialized  from  time  zero. 

a 

The  velocity  components  on  the  surface  of  the  cylinder1  were  brought ■ 
linearly  from  their  freestream  value  to  zero  in  a  period  of  6  usee 
(t^  =  6x10  ^  sec).  The  boundary  conditions  prescribed  for  this  Case 
are  those  used  with  the  finite  outer  boundary. 

This  case  was  run  from  time  zero  to  a  non-dimensional  time;  of 
t  -  4.80  at  which  time  the  flow  field  was  still  unsteady.  These 
calculations  required  approximately  40  hours  of  computer  time  on . 

BRLESC  II.  On  the  basis  of  the  shock  speed  and  standoff  distance  at 

I 

t  =  4.80  was  estimated  that  an  additional  30  to  40  hohrs  of  computer! 
time  would  be  required  before  the  flow  field  would  reach  steady-state 
conditions.  The  calculations  for  this  case  were  not  carried  beyond 
this  point. 

A  particular  difficulty  arose  with  the  density  calculations  at  the 
forward  stagnation  point  on  the  cylinder.  The  density  at  the  forward 
stagnation  point  increased  rapidly  with  time  from  its  freestream  value 
at  time  zero,  reaching  a  peak  value  approximately  12  times  the  free 4 
stream  density  at  t  =  6  usee.  Figure  4  shows  the  density  variation  in 
the  vicinity  of  the  forward  stagnation  point  for  the  first  12  usee. 

For  times  after  t  =  6  usee  the  density  at  the  forward  stagnation  point 

j 

decreased  monotonically  with  time  .  At  t  =  12  usee  the  stagnation  point 
density  had  decreased  to  approximately  seven  times  the  freestream 
value . 

A  second  calculation  was  run  for  this  case  with  t  -  3  usee,  one- 
half  the  value  used  previously  ,  The  peak  stagnation  point  density 


*  The  calculationt  fan  tklt  cate  wcAc  pen  famed  by  Vn,  Paxil  Gondon 
w fille  a  tummci  employee  at  Ike  BRL.  He  adapted  kit  original  compute a 

pnognrn^  to  BRLESC  II  and  checked  It  fan  compatibility  with  dun  computen 
by  nunnlng  a  pnevloutly  computed  cate. 
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Figure  k.  Density  Along  Stagnation  Streamline  (9  =  0C )  for 


was  lowered  by  approximately  two  percent,  but  now  occurred  at  t  =  3  nsec 
instead  of  t  =  6  (xsec.  These  two  calculations  show  that  the  time 
history  of  the  stagnation  point  density  is  quite  dependent  on  the 
initialization  process  prescribed. 

The  peak  values  of  the  density  predicted  at  the  foward  stagnation 
point  are  incorrect.  Figure  4  shows  that  there  is  a  density 
discontinuity  for  t  =••  12  |J.sec  at  the  forward,  stagnation  point.  This 
discontinuity  results  because  the  numerical  computations  do  not 
satisfy  the  symmetry  condition 

=  0  (IT) 

along  0  *  0°  .  By  additionally  specifying  (17)  this  discontinuity  is 
eliminated. 

B.  Re.  =  46.8 

d _ 

Three  sets  of  steady-state  numerical  solutions  were  obtained  for 
Re^  ~  46.8.  These  calculations  were  for  a  cylinder  of  approximately 
l/32  inch  diameter.  The  first  case  (47 -l)  was  started  at  time  zero 
and  run  until  the  steady-state  solution  was  reached  at  t  =  26.2 
after  approximately  3100  time  steps.  This  case  required  13  hours  of 
computer  time  on  BRLESC  II.  Figure  3  shows  the  placement  of  the  outer 
boundary  used  for  case  47-1. 

A  second  calculation  (47-2)  was  run  to  evaluate  the  effect  of 
specifying  freestream  conditions  at  the  outer  boundary.  For  this 
second  case  the  outer  boundary  was  placed  farther  from  the  cylinder 
than  in  the  first  case  as  is  shown  in  Figure  3-  This  case  was 
initialized  at  t  =  26.2  by  interpolating  from  the  steady-state  flow 
field  of  case  47-1.  The  governing  equations  were  then  integrated  in 
time  until  a  new  steady-state  solution  was  reached  for  this  boundary 
location.  This  required  approximately  2200  additional  time  steps. 

The  numerical  solution  achieved  a  new  steady  state  at  T  -  42.1  after 
18  additional  hours  of  computer  time. 


de 
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The  first  two  cases  (47-1  and  47-2)  used  the  "boundary  conditions 
described  previously  for  use  with  a  finite  outer  boundary.  No 
additional  density  boundary  conditions  were  specified.  The  calculations 
for  case  47-1  were  Initialized  with  t^  =  2  nsec.  A  value  of  was  not 
necessary  for  case  47-2  since  it  was  started  by  interpolation  at 
t  =26.2. 

The  third  case  (47-3)  employed  the  coordinate  mapping  (13)  with 
the  outer  boundary  at  infinity.  In  this  case  we  additionally  specified 
a  symmetry  condition  for  the  density. 


(18) 


along  the  boundaries  at  9  =  o°  and  0  =  180° .  This  case  was  started 
at  time  zero  with  t  <  At  and  thus  treated  a  cylinder  impulsively 
accelerated  from  rest.  The  steady-state  solution  was  reached  at 
t  =  39*98  after  12  hours  of  computer  time  and  approximately  2900 
time  steps. 

None  of  the  steady-state  solutions  for  Re^  =  46. S  predicted 
separation.  The  effect  of  the  downstream  symmetry  conditions  on  the 
numerical  solutions  could  not  be  determined  from  these  calculations. 
In  order  to  determine  if  the  symmetry  boundary  conditions  were 
restricting  "numerical  separation"  a  fourth  case  47-4  was  run. 

Case  47-4  was  set-up  with  a  difference  mesh  surrounding  the 
cylinder  (0  ^  0  £  360°).  The  number  of  mesh  elements  used  was  twice 
that  of  case  47-3.  The  mesh  size  was  identical  to  that  used  in  case 
47-3.  This  360°  placement  of  nodes  eliminated  the  necessity  for  the 
specification  of  downstream  symmetry  conditions.  This  case  was 
started  at  time  zero  in  the  same  manner  as  case  47-3  and  run  to  a 
non-dimensional  time  t  =  7.90.  A  comparison  of  cases  47-3  and  47-4 
showed  them  to  be  almost  identical  at  each  time  step  and  it  was 
concluded  that  their  steady-state  solutions  would  be  the  same.  The 
calculations  for  case  47-4  were  stopped  at  this  point.  These  results 
indicated  that  separation  was  not  restricted  by  our  symmetry  boundary 
conditions. 


28 


Figure  5  shows  a  comparison  of  the  three  steady-state  numerical 
solutions  in  the  wake  region.  Near  the  cylinder  (r/r  <  3)  the  two 
calculations  using  a  finite  outer  boundary  (oases  47-1  and  47-2)  are 
in  agreement.  The  calculations  for  case  47-3,  using  the  coordinate 
mapping i  predict  higher  densities  near  the  cylinder  than  predicted, 
by  the  first  two  cases.  The  calculations  for  cases  47-1  and  47-2 
fail  far  downstream  of  the  cylinder.  These  two  cases  predict  large 
gradients  of  velocity,  density  and  temperature  far  downstream.  In 
the  physical  flow  field  large  gradients  would  be  present  if  there  were 
a  wake  shock  but  would  be  opposite  in  sign  to  those  predicted  by  these 
calculations.  The  large  gradients  predicted  downstream  in  cases  47-1 
and  47-2  are  typical,  of  those  in  rarefaction  waves.  These  gradients 
result  from  the  unnatural  acceleration  the  flow  experiences  in  the 
numerical  solution  by  the  specification  of  freestream  boundary 
conditions  downstream  of  the  cylinder. 

The  numerical  solution  for  case  47-3  predicts  a  more  realistic 
qualitative  behavior  far  downstream  (r/Rq  >  10).  The  large  gradients 
are  eliminated  in  this  solution  as  a  result  of  the  asymptotic  behavior 
of  the  coordinate  mapping .  None  of  the  three  steady-state  solutions 
for  Re^  =  46.8  predicted  a  wake  shock  downstream  of  the  cylinder. 

In  cases  47-1  and  47-2  no  density  boundary  condition  was  specified 
along  the  downstream  portion  of  the  outer  'boundary  (refer  to  4  in 
Figure  3).  Both  of  these  cases  predict  densities  much  lower  than  free¬ 
stream  at  the  outer  boundary  along  0  =  l80° .  Case  47-3  predicts  a 
density  downstream  which  asymptotically  approaches  the  freestream 
value.  It  thus  appears  that  the  mapping  and  boundary  conditions  used 
in  case  ^7-3  produce  better  qualitative  results  downstream,  of  the 
cylinder  than  the  calculations  prescribing  freestream  boundary 
conditions  along  the  finite  outer  boundary. 

Figure  6  shows  the  density  distribution  across  the  bow  shock,  at 
its  intersection  with  the  stagnation  streamline.  The  numerical 
solutions  for  Re^  -  46.8  predict  shock  thicknesses  of  two  to  three 
mean  free  paths.  The  analytical  viscous  normal  shock  solution  of 
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Morduchow  and  Libby‘S  shown  in  this  figure  has  a  shock  thickness  of 
two  mean  free  paths.  The  good  resolution  of  the  shock  region  in  our 
numerical  computations  is  due  to  the  very  fine  difference  mesh  used 
in  the  bow  shock  region j  less  than  one  mean  free  path  in  the  radial 
direction.  Numerical  solutions  typically  spread  shocks  over  many 

12 

mesh  widths  and  thus  predict  innacurate  shock  thicknesses.  Walitt 
calculates  shock  thicknesses  on  the  order  of  90  mean  free  paths  for 
Mach  two  flow  at  Re^  =  l4Q9* 

Figure  7  shows  the  time  history  of  the  density  ratio  across  the 

bow  shock  at  9  =  0°  .  The  numerical  calculations' for  case  47-1  predict 

a  steady-state  density  ratio  which  is  seven  percent  higher  than  that 

predicted  analytically.  The  analytical  solution  used  for  comparison 

is  the  normal  shock  solution  with  first-order  corrections  for  shock 

13 

curvature  due  to  Chow  and  Ting  .  The  numerical  calculations  for  case 

47-3  predict  a  steady-state  density  ratio  which  is  14  percent  higher 

than  this  analytical  solution.  The  Rankine -Hugoniot  normal  shock 

solution  is  also  shown  in  Figure  7. 

Figure  8  shows  the  density  distribution  along  the  stagnation 

streamline  predicted  by  the  three  steady-state  numerical  solutions. 

Near  the  cylinder  case  47-3  predicts  much  higher  densities  than 

either  case  47-1  or  47-2.  Case  47-2  is  seen  to  predict  a  larger 

standoff  distance  than  either  of  the  other  two  cases.  These  results 

show  that  the  manner  in  which  the  boundary  conditions  are  specified 

for  the  numerical  solutions  has  a  very  large  influence  on  the  density 

distribution  predicted  near  the  cylinder j  both  upstream  and  downstream. 

None  of  the  three  steady-state  calculations  for  Re^  =  46.8 

showed  any  tendency  toward  separation  as  we  would  expect  on  the  basis 

14 

of  the  incompressible  experimental  results  by  Taneda  .  His  results 
show  that  for  incompressible  flow  past  cylinders  a  separation  bubble 
forms  downstream  for  Reynolds  numbers  above  Re^  «  5.  For  Reynolds 
numbers  above  Re^  **  45  part  of  the  separation  bubble  is  shed  alternately 
from  each  of  the  two  recirculation  regions.  In  the  physical 
incompressible  flow  field  symmetric  separation  is  unstable  for 


32 


Density  Distribution  Along  Stagnation  Streamline 


Reynolds  numbers  above  Re^  sa  4-5.  A  calculation  for  Re^  -15.0  was 
run  to  make  certain  that  our  symmetric  downstream  boundary  conditions 
were  compatible  with  the  stability  of  the  flow  field. 

C.  Red  -  15.0 

The  numerical  calculations  for  Re^  =  15*0  employed  the  coordinate 
mapping  (13)  and  the  boundary  conditions  previously  described  for  use 
with  the  mapping.  The  computations  were  started  at  time  zero  using 
the  initial  conditions  for  the  impulsively  accelerated  cylinder.  The 
steady-state  solution  was  reached  at  t  =47.7  after  13  hours  of 
computer  time  and  approximately  3000  time  steps. 

The  mesh  size  used  in  this  computation  was  identical  to  that 
used  for  case  47-3*  The  density  was  the  only  freestream  flow 
variable  which  differed  from  that  used  in  case  47-3 •  Table. AI  shows 
the  values  specified  for  the  freestream  flow  quantities  in  this  case. 

The  purpose  of  this  low  Reynolds  number  ease  was  to  determine  if 
our  specification  of  downstream  symmetry  boundary  conditions  restricted 
separation  in  the  higher  Reynolds  number  solutions  described  previously. 
The  steady-state  numerical  solution  for  case  15  showed  no  tendency 
toward  separation.  It  thus  appears  that  our  symmetry  boundary  condition 
is  compatible  with  the  numerical  computation  and  does  not  prevent 
separation. 

Figure  9  shows  the  steady-state  density  and  temperature  distri¬ 
butions  along  the  stagnation  streamline  for  Rc^  =  15*0.  The  density 
ratio  predicted  by  this  numerical  solution  across  the  bow  shock  is 

2.77.  This  is  three  percent  higher  than  that  predicted  by  the  normal 

13 

shock  solution  with  first-order  corrections  for  shock  curvatuve 
The  standoff  distance  predicted  by  this  numerical  solution  is 
expectantly  larger,  approximately  45  percent  larger  than  that 
predicted  by  the  numerical  calculations  for  Re^  =  46.8. 

It  is  interesting  to  compare  this  solution  for  case  15  with  the 
solution  for  case  47-3.  Both  of  these  cases  employed  the  same  mesh 
sizes  with  all  freestream  flow  quantities  identical  except  for  the 
density.  Cases  15  and  47-3  used  identical  initial  conditions  and 
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boundary  conditions  and  both  used  the  coordinate  mapping  (15) •  The 
numerical  solutions  in  the  bow  shock  region  are  quite  different. 

Case  47-3  predicts  a  density  ratio  across  the  bow  shock  which  is 
much  too  high  as  described  previously.  Case  15  predicts  a  density 
ratio  across  the  bow  shock  which  is  quite  close  to  the  analytical 
result. 

Since  both  case  15  and  case  47-5  use  the  same  initial  conditions , 
mapping,  and  boundary  conditions  it  appears  that  the  number  of  mesh 
points  across  the  bow  shock  is  the  crucial  difference  between  these 
two  calculations.  Case  15  has  almost  three  times  as  many  mesh  points 
across  the  bow  shock  at  9  =  0°  as  case  47-5*  This  is  because  the  gas 
mean  free  path  is  larger  in  case  15  and  the  shock  is  consequently 
thicker.  Case  15  had  12  mesh  points  across  the  bow  shock  where  it 
intersected  the  stagnation  streamline  (almost  four  mesh  points  per 
mean  free  path) . 

D.  Red  =  157-2 

A  numerical  calculation  was  made  for  a  Reynolds  number  of  Re^  = 
157.2.  The  purpose  of  this  calculation  was  to  obtain  a  numerical 
solution  which  could  be  compared  with  experimental  wake  measurements 
taken  for  the  same  Mach  number  and  Reynolds  number.  These  hot-wire 
measurements  taken  downstream  of  a  circular  cylinder  at  Re^  =  157.2 
in  a  Mach  two  freestream  show  that  a  wake  shock  is  present  in  the 
flow  field. 

Several  attempts  were  made  to  start  the  integration  of  case  157 
at  time  zero  using  the  "impulsive  acceleration'1  initial  conditions. 
These  initial  conditions  produced  excessive  density  ratios  at  the 
forward  stagnation  point  on  the  cylinder.  The  density  ratio  reached 
values  in  excess  of  an  order  of  magnitude  larger  than  that  for  a 
flow  undergoing  an  unsteady  isentropic  compression.  Due  to  these 
inaccurate  numerical  results  during  the  first  few  time  steps  this 
procedure  was  abandoned  in  favor  of  interpolation  from  another  steady- 
state  solution. 
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The  numerical  computation  of  case  157  was  Btarted  at  t  =  39.08  by 
Interpolating  from  the  steady-state  flow  field  of  case  1+7-3.  The  free- 
stream  density,  velocity,  and  temperature  used  for  case  157  were 
identical  to  that  used  in  case  1+7-3.  The  viscosity  relation  used  for 
case  157  was  varied  from  that  used  in  case  47-3  to  produce  a  higher 
Reynolds  number.  Table  AI  shows  the  values  used  to  prescribe  the  free- 
stream  flow  field. 

The  numerical  solution  reached  the  steady-state  at  T  =  58.14 
after  nine  hours  of  computer  time  and  approximately  1500  time  steps. 
Since  the  initial  computations  were  performed  by  interpolating  from 
a  lower  Reynolds  number  solution  the  time  dependence  for  case  157  does 
not  represent  a  cylinder  accelerated  from  rest.  The  use  of  this 
interpolation  technique  for  starting  a  solution  offers  a  possible 
approach  for  obtaining  an  even  higher  Reynolds  number  solution  in  a 
reasonable  amount  of  computer  time. 

The  numerical  computation  for  Re^  =  157-2  predicted  a  bow  shock 
standoff  distance  of  A/Rq  =  1.282.  Figure  10  shows  a  comparison  of  the 
time  history  of  the  standoff  distance  for  each  calculation.  The  stand¬ 
off  distance  calculated  for  Re,  =  157.2  is  six  percent  larger  than  that 

d 

predicted  by  inviscid  calculations.  Our  numerical  solutions  predict 
that  the  3tandoff  distance  increases  with  decreasing  Reynolds  number. 
This  variation  of  standoff  distance  is  as  we  would  expect. 

The  density  ratio  across  the  bow  shock  predicted  by  the  numerical 
calculations  for  Re^  =  157-2  is  much  too  large.  At  the  point  where 
the  how  shock  intersects  the  stagnation  streamline  a  density  ratio 
of  p0/p  =  3-5  is  predicted  by  the  numerical  computations.  This  is 
29  percent  larger  than  predicted  analytically  by  viscous  normal  shock 
solutions  corrected  to  account  for  shock  curvature"^  The  cause  of 
this  large  discrepancy  appears  to  be  that  more  mesh  points  than  we 
used  are  needed  in  the  vicinity  of  the  bow  shock. 

Figure  11  shows  the  density  distribution  through  the  bow  shock 
along  8=0°  in  terms  of  the  freestream  mean  free  path.  The  numerical 
calculations  for  Re^  =  157-2  predict  a  shock  thickness  of  approximately 
three  mean  free  paths  while  the  analytical  viscous  normal  shock 
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Figure  10.  'Time  History  of  Bow  Shock  Standoff  Distance 
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solution  predicts  a  shock  thickness  of  two  mean  free  paths.  Also 
apparent  in  Figure  11  is  the  large  discrepancy  between  the  Rankine- 
Hugonlot  and  numerical  solution  for  the  density  ratio  across  the  bow 
shock . 

Figure  12  shows  the  pressure  coefficient  predicted  throughout  the 

flow  field  for  case  157-  Upstream  of  the  cylinder  (  X/Rq  =  -  2)  the 

bow  shock  is  quite  pronounced  By  the  time  the  flow  reaches  X/Rq  =  1 

the  bow  shock  has  been  smeared  over  many  mesh  widths.  This  smearing 

or  weakening  of  the  bow  shock  is  seen  to  continue  as  the  flow  moves 

downstream.  The  pressure  coefficient  at  the  forward  stagnation  point 

is  C  =  2.29.  The  bow  shock  is  being  smeared  out  by  the  asymptotic 
1? 

behavior  of  the  coordinate  mapping  (13 )  used  for  the  computations. 

This  smearing  or  weakening  of  the  bow  shock  is  caused  because  free- 
s bream  conditions  are  specified  at  infinity  and  the  condition 

d/3R  “*OasR-*<D  (l?) 

is  inherent  to  the  coordinate  mapping  (13).  The  approximate  location 
of  the  wake  shock  is  also  shown  in  Figure  12  for  reference. 

Figure  13  shows  the  pressure  coefficient  in  the  wake  region 
predicted  by  the  numerical  solution  for  case  157-  Near  the  rear 
stagnation  point  (0  =  l80° )  the  pressure  is  less  than  its  freestream 
level.  By  the  time  the  flow  along  the  rear  stagnation  streamline 
reaches  X/Rq  ~  4  its  pressure  level  has  exceeded  that  in  the  freestream. 
These  profiles  exhibit  "shock-like"  gradients  in  the  wake  region 
(i.e.,  X/Rq  =  3,  Y/Rq  ~  1.5).  The  dashed  curve  shown  in  Figure  13 
represents  the  locus  of  points  of  maximum  positive  pressure  gradient. 
Pressure  increases  as  we  cross  this  dashed  line  moving  downstream. 

These  gradients  result  from  a  compression  of  the  flow  and  appear  to 
represent  the  woke  shock  in  the  physical  flow  field.  The  numerical 
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and  experimental  results  in  the  wake  region  cannot  be  compared 
directly  since  the  experimental  measurements  were  taken  at  x/Rq  70 
and  the  wake  shock  is  smeared  out  in  our  numerical  solution  by 
X/Rq  ~  13  due  to  our  placement  of  node  points  and  the  asymptotic 
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behavior  of  the  coordinate  mapping.  We  conclude  that  the  numerical 
computations  predict  "shock -like 11  gradients  in  the  wake  at  Re^  =  157.2 
which  may  represent  the  wake  shock  present  in  the  physical  flow  field 
at  this  same  Reynolds  number. 

The  numerical  computations  for  Re^  =  157.2  predicted  a  small 
recirculation  region  downstream  of  the  cylinder.  The  separation 
bubble  predicted  in  these  numerical  calculations  is  shown  in  Figure  14. 
Separation  is  predicted  to  start  at  6  ~  1J80.  On  the  basis  of 
incompressible  flow  results  this  separation  is  very  much  delayed.  For 
incompressible  flow  the  separation  point  is  located  at  0  — •  109° •  At 
the  present  time  it  is  hot  known  whether  this  delayed  separation  is  a 
result  of  compressibility  effects  or  errors  in  the  numerical 
computations.  It  should  be  noted  tliat  the  numerical  computations  for 
Re^  -  15.0  and  46.8  predicted  no  separation  bubble  downstream  of  the 
cylinder. 

VI.  CONCLUSIONS 

This  report  lias  described  a  series  of  numerical  experiments  for 
the  purpose  of  evaluating  the  accuracy  of  typical  numerical  solutions 
of  the  full  Navier- Stokes  equations.  A  computer  program  developed  by 
Scala  and  Gordon'1  was  used  in  this  study  to  solve  for  the  steady- 
state  flow  field  around  cylinders  accelerated  to  supersonic  velocities. 
Steady-state  solutions  were  obtained  for  Mach  two  flow  past  cylinders 
at  Reynolds  numbers  of  15.O,  46.8,  and  157.2. 

Oui’  computations  have  shown  that  the  density  levels  predicted 

downstream  of  the  bow  shock  are  higher  than  we  would  expect  on  the 

basis  of  analytical  normal  shock  solutions  corrected  to  account  for 

shock  curvature1'^.  The  degree  of  error  present  in  our  numerical 

solutions  is  closely  linked  to  the  spacing  of  mesh  points  in  regions 

where  large  gradients  are  present.  The  steady-state  calculations  for 

Re,  =  15.0  had  12  mesh  points  across  the  bow  shock  and  had  three  per- 
d 

cent  toe  high  a  density  ratio;  the  calculation  for  Re^  =  157*2  had 
29  percent  tco  high  a  density  ratio  but  only  5  mesh  points  across 
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14.  Separation  Bubble  in  numerical  Calculations 


the  bow  shock.  It  thus  appears  necessary  to  UBe  very  fine  mesh 

opacings  in  flow  regions  having  large  gradients  for  accurate 

numerical  calculations.  Three  or  four  mesh  points  per  mean  free 

path  appear  to  be  necessary  in  the  region  of  the  bow  shock. 

The  numerical  solutions  for  Re^  =  15*0  and-  ^6.8  did  not 

predict  separation.  The  computation  for  Re^  =  157.2  predicted  a 

small  separation  bubble  on  the  downstream  side  of  the  cylinder.  For 

this  case  separation  began  at  0  =  180°  then  the  separation  point 

shifted  with  time  reaching  a  steady-state  position  at  9~  158°  Our 

numerical  results  thus  predict  very  late  separation  compared  to 

incompressible  flow  past  a  cylinder .  Walitt  predicts  separation  at 

0  ~  i4l°  in  a  numerical  calculation  for  Mach  two  flow  past  a  cylinder 

at  Re  =  1409.  The  accuracy  of  our  numerically  predicted  separation 
d 

point  is  unknown. 

Steady-state  numerical  solutions  were  obtained  for  Re^  =  15*0> 

46.8  and  157.2  by  employing  the  coordinate  mapping  suggested  by 
Sills'10.  Our  calculations  have  shown  that  this  coordinate  mapping  is 
compatible  with  initial -boundary-value  problems  of  the  present  type. 
This  mapping  smooths  out  gradients  far  from  the  body.  This  smoothing 
process  is  beneficial  in  the  wake  region.  The  mapping  eliminated  the 
large  unnatural  gradients  which  resulted  from  specifying  artificial 
freoatream  conditions  downstream  of  the  cylinder.  The  two  methods  of 
specifying  boundary  conditions  downstream  are  discussed  and  compared 

in  the  calculations  for  Rc^  •-  46.8. 

Experimental  wake  measurements  at  Re^  =  157-2  have  shown  that  a 

wake  shock  is  present  in  Mach  two  flow  at  this  Reynolds  number.  Our 

numerical  computations  for  Re^  =  157*2  predict  large  pressure  and 

temperature  gradients  at  the  edge  of  the  inner  wake  which  appear  to  be 

the  numerical  equivalent  of  the  wake  shock  in  the  physical  flow  field. 

It  is  not  clear,  however,  that  the  details  of  the  wake  region  are 

accurately  predicted  by  these  calculations.  The  numerical  calculations 

for  Re  =  15.0  and  46.8  predicted  no  wake  shock.  There  is  no  experi- 
d 

mental  data  for  Reynolds  numbers  as  low  as  46.8. 


The  present  numerical  experiments  have  been  for  cylinders  of 
approximately  l/j2  inch  and  l/2  inch  diameter.  To  maintain  a  desired 
level  of  accuracy  the  number  of  mesh  points  used  in  a  computation  must 
increase  if  larger  cylinders  are  considered.  In  view  of  the  large 
computer  times  required  for  our  calculations  these  small  cylinder 
sizes  represent  an  upper  limit  of  our  present  capability.  The 
calculation  of  a  typical  steady-state  solution  required  from  nine  to 
eighteen  hours  on  the  BRLESC  II  computer.  For  a  given  problem  it  is 
impossible  to  predict  in  advance  the  degree  of  error  which  will  be 
present  in  the  steady-state  numerical  solution.  Our  computations  have 
shown  that  errors  of  ten  to  twenty  percent  may  result  in  numerical 
solutions  of  the  Navier -Stokes  equations  even  under  carefully 
controlled  conditions. 
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APPENDIX  A 


The  following  gas  properties  were  held  fixed  for  all  the 
numerical  cases  considered 

R  =  5-52X101  fft-lbj/lfc  -  °R  ,  (Al) 

a  =  1.38x102  ft  -lb./lb  -  °R  ,  (AS) 

v  r  m 

Y  =  1.40  .  (A3) 

The  freestream  Reynolds  number  was  varied  by  changing  the  cylinder 
size  or  freestream  conditions  of  the  flow.  The  specific  values  used 
for  each  case  are  shown  in  Table  A -I.  These  dimensional  variables 
are  input  quantities  for  the  computer  program  developed  by  Scala  and 
Gordon"*".  Each  case  is  specified  by  its  approximate  Reynolds  number. 
The  four  cases  for  Re^  =  46.8  are  differentiated  by  a  numerical 
suffix  (i.e.  47-l).  Table  A -I  also  shows  which  calculations  were 
performed  using  the  coordinate  mapping  and  the  number  of  meBh  points 
used  for  each  calculation. 


Preceding  page  blank 
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Dimensional  Values  of  Free stream  Variables  and  Size  of  Computation 


